20 October 2017

Why?

Why use R as a GIS?

1- Workflow

2- Quite efficient

3- Interface

Workflow

R for everything (almost)

  • import your data
  • format your data
  • analyze your data
  • visualize your data
  • export your data
  • create your own function/package

Efficiency

  • well-defined classes
  • many formats (+ convert different file)
  • topology operations
  • statistical analyses
  • visualization

Interface

Overview

Turning R into a powerful GIS

1- Define classes

  • sp
  • raster

2- Topology operations

  • rgeos

3- Import / export

  • rgdal

Turning R into a powerful GIS

5- Analyses

  • spdep
  • dismo

4- Creating and editing maps

  • graphics
  • ggplot
  • interactive plot

6- But

  • Geo-referencing …

Package sp

Classes

Classes / Functions Contents
Points list of points (set of coordinates)
SpatialPoints list of points + CRS
SpatialPointsDataPoints list of points + CRS + attribute table
Line a line (set of coordinates)
Lines list of lines
SpatialLines list of lines + CRS
SpatialLinesDataFrame list of lines + CRS + attribute table

Example: SpatialPointsDataPoints

long lat var1 var2
-80.55252 43.51967 -0.3811779 0.7639132
-81.85050 43.25344 -1.3384095 9.2054164
-80.60760 42.14834 0.6091255 5.3426150
-80.23205 42.70182 -0.6528624 3.1712027
-81.47843 42.11437 0.5111311 3.3774720
-80.88094 42.61269 2.0061587 7.7422934

Example: SpatialPointsDataFrame

library(sp)
## Loading required package: methods
mysp <- SpatialPointsDataFrame(
  coords = mydata[,1:2],
  data = mydata[,3:4],
  proj4string = CRS(
    "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"
  )
)

Example: SpatialPointsDataFrame

class(mysp)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
isS4(mysp)
## [1] TRUE
slotNames(mysp)
## [1] "data"        "coords.nrs"  "coords"      "bbox"        "proj4string"

Example: SpatialPointsDataFrame

mysp@proj4string
## CRS arguments:
##  +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0
head(mysp@data)
##         var1      var2
## 1 -0.3811779 0.7639132
## 2 -1.3384095 9.2054164
## 3  0.6091255 5.3426150
## 4 -0.6528624 3.1712027
## 5  0.5111311 3.3774720
## 6  2.0061587 7.7422934

Change projection: spTransform

(mysp2 <- spTransform(mysp, CRS=CRS("+proj=merc +ellps=GRS80")))
##            coordinates        var1      var2
## 1  (-8967066, 5361976) -0.38117786 0.7639132
## 2  (-9111556, 5321339) -1.33840953 9.2054164
## 3  (-8973197, 5154544)  0.60912553 5.3426150
## 4  (-8931391, 5237708) -0.65286244 3.1712027
## 5  (-9070137, 5149464)  0.51113113 3.3774720
## 6  (-9003625, 5224266)  2.00615870 7.7422934
## 7  (-9024502, 5259896) -2.16004226 9.1975772
## 8  (-8973962, 5320027)  0.24658562 1.3054099
## 9  (-9029957, 5344375) -0.17591600 4.5683469
## 10 (-9121643, 5324546) -1.58085174 5.4410930
## 11 (-9086762, 5185078) -0.95839638 1.5291094
## 12 (-9003765, 5415044) -1.02981222 4.7838512
## 13 (-9059708, 5307232)  0.52491399 6.7900775
## 14 (-8982200, 5161679)  0.76628680 4.9152753
## 15 (-9051522, 5314610)  0.21809976 3.4204696
## 16 (-9029497, 5294779) -1.09869789 3.4025504
## 17 (-9051360, 5160753)  0.01876083 2.3815339
## 18 (-9073265, 5308913)  0.05346164 7.8725486
## 19 (-8946979, 5156041)  1.39310716 3.9224266
## 20 (-9103984, 5208464) -1.11292896 3.7748791

Package raster

Classes

1- SpatialGrid (package sp)

2- SpatialGrid (package sp)

3- RasterLayer (package raster)

4- RasterStack (package raster)

5- RasterBrick (package raster)

Example:raster

library(raster)
## 
## Attaching package: 'raster'
## The following object is masked from 'package:magrittr':
## 
##     extract
ras1 <- raster(matrix(runif(100*100,0,10),ncol=100,nrow=100),
    crs=CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"),
    xmn=-82, xmx=+80, ymn=+42, ymx=+44)

Example:raster

ras1
## class       : RasterLayer 
## dimensions  : 100, 100, 10000  (nrow, ncol, ncell)
## resolution  : 1.62, 0.02  (x, y)
## extent      : -82, 80, 42, 44  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : layer 
## values      : 0.00157753, 9.999279  (min, max)

Retrieving free GIS data: getData

head(getData("ISO3"))
##   ISO3        NAME
## 1  ABW       Aruba
## 2  AFG Afghanistan
## 3  AGO      Angola
## 4  AIA    Anguilla
## 5  ALA       Åland
## 6  ALB     Albania

Example:getData

## Country level:
mapBEL0 <- getData(name="GADM", country="BEL", path="./assets", level=0)
mapBEL1 <- getData(name="GADM", country="BEL", path="./assets", level=1)
tminW <- getData(name="worldclim", var="tmin", res=10, path="./assets")
mapBEL0
## class       : SpatialPolygonsDataFrame 
## features    : 1 
## extent      : 2.555356, 6.40787, 49.49722, 51.50382  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 68
## names       : OBJECTID, ID_0, ISO, NAME_ENGLISH, NAME_ISO, NAME_FAO,              NAME_LOCAL, NAME_OBSOLETE,            NAME_VARIANTS, NAME_NONLATIN, NAME_FRENCH, NAME_SPANISH, NAME_RUSSIAN, NAME_ARABIC, NAME_CHINESE, ... 
## min values  :        1,   23, BEL,      Belgium,  BELGIUM,  Belgium, Belgique|Belgie|Belgien,              , Belgique|Belgien|Belgium,              ,   Belgique ,     Bélgica ,     Бельгия ,      بلجيكا,      比利时 , ... 
## max values  :        1,   23, BEL,      Belgium,  BELGIUM,  Belgium, Belgique|Belgie|Belgien,              , Belgique|Belgien|Belgium,              ,   Belgique ,     Bélgica ,     Бельгия ,      بلجيكا,      比利时 , ...

Package rgdal

Drivers

library(rgdal)
## rgdal: version: 1.2-13, (SVN revision 686)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.1.2, released 2016/10/24
##  Path to GDAL shared files: /usr/share/gdal/2.1
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: (autodetected)
##  Linking to sp version: 1.2-5
  1. writeOGR / writeGDAL: to write spatial objects
  2. readOGR/readGDAL: to read spatial files

Drivers

head(ogrDrivers())
##         name                     long_name write  copy isVector
## 1 AeronavFAA                   Aeronav FAA FALSE FALSE     TRUE
## 2 AmigoCloud                    AmigoCloud  TRUE FALSE     TRUE
## 3     ARCGEN             Arc/Info Generate FALSE FALSE     TRUE
## 4     AVCBin      Arc/Info Binary Coverage FALSE FALSE     TRUE
## 5     AVCE00 Arc/Info E00 (ASCII) Coverage FALSE FALSE     TRUE
## 6        BNA                     Atlas BNA  TRUE FALSE     TRUE

Export

writeOGR(mysp, dsn="./assets", layer="mypoints",
    driver="ESRI Shapefile", overwrite_layer=TRUE)

Import

mysp2 <- readOGR(dsn="assets/", layer="mypoints")
## OGR data source with driver: ESRI Shapefile 
## Source: "assets/", layer: "mypoints"
## with 20 features
## It has 2 fields
## Roads find on http://www.diva-gis.org/Data
canroads <- readOGR(dsn="assets/", layer="CAN_roads")
## OGR data source with driver: ESRI Shapefile 
## Source: "assets/", layer: "CAN_roads"
## with 16915 features
## It has 5 fields

Note: there are also functions to import/export raster in raster

Package rgeos

Load the package

library(rgeos)
## rgeos version: 0.3-25, (SVN revision 555)
##  GEOS runtime version: 3.5.1-CAPI-1.9.1 r4246 
##  Linking to sp version: 1.2-5 
##  Polygon checking: TRUE

Load the package

buf <- gBuffer(mapBEL0, width=0.5)
## Warning in gBuffer(mapBEL0, width = 0.5): Spatial object is not projected;
## GEOS expects planar coordinates
diff <- gDifference(gBuffer(mapBEL0, width=0.5), gBuffer(mapBEL0, width=0.1))
## Warning in gBuffer(mapBEL0, width = 0.5): Spatial object is not projected;
## GEOS expects planar coordinates
## Warning in gBuffer(mapBEL0, width = 0.1): Spatial object is not projected;
## GEOS expects planar coordinates

Package mapview

Import package

library(mapview)
## Loading required package: leaflet

NB: it uses the leaflet package.

Quick examples

mapview(mysp, cex = 'var2')@map

Quick examples

mapview(mapBEL1)@map

More

Editing a map

A very basic map – Shapefile

plot(mapBEL0)

A very basic map – Raster

class(tminW)
## [1] "RasterStack"
## attr(,"package")
## [1] "raster"

A very basic map – Raster

plot(tminW)

Cutomize a map – Shapefile

plot(mapBEL0, border='grey15', col='#E6E6E6', lwd=1.6)
plot(mapBEL1, lty=2, lwd=0.9, add=T)
points(4.3513, 50.8471, pch=19, col="#27df9d")
text(4.3513, 50.8471, text="Brussel", pos=3)

Cutomize a map – Shapefile

Cutomize a map – Shapefile of roads

Examples

Data frame

  • this .csv is avaible on line at XXXX
lakedf <- read.csv('assets/lakeOnt.csv')
head(lakedf)
##   X LAKEID LONGITUDE_DD LATITUDE_DD Area_ha
## 1 1 L32173    -80.03361    49.90417   174.9
## 2 2 L21001    -91.68278    48.68028   130.2
## 3 3 L44001    -84.19667    46.91861   279.6
## 4 4 L25002    -90.43111    48.16639   114.1
## 5 5 L44002    -84.30528    47.06000   137.3
## 6 6 L21003    -91.34528    48.24806  2993.9

SpatialPointDataFrame

library(sp)
lakesp <- SpatialPointsDataFrame(
  coords = lakedf[,3:4],
  data = lakedf[,c(1,4)],
  proj4string = CRS(
    "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"
  )
)

SpatialPointDataFrame

lakesp
## class       : SpatialPointsDataFrame 
## features    : 546 
## extent      : -95.06111, -76.01194, 43.05722, 54.52056  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0 
## variables   : 2
## names       :   X, LATITUDE_DD 
## min values  :   1,    43.05722 
## max values  : 546,    54.52056

Altitude Canada

altCAN <- getData(name="alt", country="CAN", path="./assets/")
altCAN
## class       : RasterLayer 
## dimensions  : 4992, 10620, 53015040  (nrow, ncol, ncell)
## resolution  : 0.008333333, 0.008333333  (x, y)
## extent      : -141.1, -52.6, 41.6, 83.2  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## data source : /home/kevcaz/Github/Tutorials/mapsWithR/docs/assets/CAN_msk_alt.grd 
## names       : CAN_msk_alt 
## values      : -108, 5800  (min, max)

Altitude Canada

Canada – provinces

bouCAN <- getData(country='CAN', level=1, path="./assets/")
bouCAN
## class       : SpatialPolygonsDataFrame 
## features    : 15 
## extent      : -141.0069, -52.61889, 41.67693, 83.11042  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 13
## names       : OBJECTID, ID_0, ISO, NAME_0, ID_1,  NAME_1, HASC_1, CCN_1, CCA_1,     TYPE_1, ENGTYPE_1, NL_NAME_1,                                       VARNAME_1 
## min values  :        1,   42, CAN, Canada,    1, Alberta,  CA.AB,    NA,    10,   Province,  Province,          ,                                                 
## max values  :       15,   42, CAN, Canada,   13,   Yukon,  CA.YT,    NA,    62, Territoire, Territory,          , Yukon Territory|Territoire du Yukon|Yukon|Yuk¢n

Canada – provinces

bouCAN@data[,1:8]
##    OBJECTID ID_0 ISO NAME_0 ID_1                    NAME_1 HASC_1 CCN_1
## 1         1   42 CAN Canada    1                   Alberta  CA.AB    NA
## 2         2   42 CAN Canada    2          British Columbia  CA.BC    NA
## 3         3   42 CAN Canada    3                  Manitoba  CA.MB    NA
## 4         4   42 CAN Canada    4             New Brunswick  CA.NB    NA
## 5         5   42 CAN Canada    5 Newfoundland and Labrador  CA.NF    NA
## 6         6   42 CAN Canada    6     Northwest Territories  CA.NT    NA
## 7         7   42 CAN Canada    7               Nova Scotia  CA.NS    NA
## 8         8   42 CAN Canada    8                   Nunavut  CA.NU    NA
## 9         9   42 CAN Canada    8                   Nunavut  CA.NU    NA
## 10       10   42 CAN Canada    8                   Nunavut  CA.NU    NA
## 11       11   42 CAN Canada    9                   Ontario  CA.ON    NA
## 12       12   42 CAN Canada   10      Prince Edward Island  CA.PE    NA
## 13       13   42 CAN Canada   11                    Québec  CA.QC    NA
## 14       14   42 CAN Canada   12              Saskatchewan  CA.SK    NA
## 15       15   42 CAN Canada   13                     Yukon  CA.YT    NA

Canada – Ontario

bouCAN[11,]
## class       : SpatialPolygonsDataFrame 
## features    : 1 
## extent      : -95.15598, -74.34605, 41.67693, 56.86972  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 13
## names       : OBJECTID, ID_0, ISO, NAME_0, ID_1,  NAME_1, HASC_1, CCN_1, CCA_1,   TYPE_1, ENGTYPE_1, NL_NAME_1,    VARNAME_1 
## min values  :       11,   42, CAN, Canada,    9, Ontario,  CA.ON,    NA,    35, Province,  Province,          , Upper Canada 
## max values  :       11,   42, CAN, Canada,    9, Ontario,  CA.ON,    NA,    35, Province,  Province,          , Upper Canada

Canada – Ontario

Canada – Ontario

Canada – Ontario elevation

altONT <- rasterize(bouCAN[11,], crop(altCAN, bouCAN[11,]@bbox), mask=TRUE)
mapview(altONT)+mapview(lakesp)

Canada – Ontario elevation

Canada – Ontario elevation

par(mar=c(0,0,0,0))
plot(altONT)
plot(lakesp, add=T)

Canada – Ontario elevation

Resources

Useful links